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In important early work, Stell showed that one can determine the pair correlation function h{r) 
of the hard sphere fluid for all distances r by specifying only the "tail" of the direct correlation 
function c(r) at separations greater than the hard core diameter. We extend this idea in a very 
natural way to potentials with a soft repulsive core of finite extent and a weaker and longer ranged 
tail. We introduce a new continuous function T{r) which reduces exactly to the tail of c(r) outside 
the (soft) core region and show that both h{r) and c(r) depend only on the "out projection" of 
r(r): i.e., the product of the Boltzmann factor of the repulsive core potential times T{r). Standard 
integral equation closures can thus be reinterpreted and assessed in terms of their predictions for the 
tail of c(r) and simple approximations for its form suggest new closures. A new and very efficient 
variational method is proposed for solving the Ornstein-Zernike equation given an approximation 
for the tail of c. Initial applications of these ideas to the Lennard- Jones and the hard core Yukawa 
fluid are discussed. 



I. INTRODUCTION 

One of the many areas of current research where George Stell has made fundamental contributions is the derivation 
of integral equations to determine the pair correlation function of a uniform fluid. A number of different integral 
equations have been proposed [1], often based on the graphical and functional methods pioneered by Stell [2]. 
However, despite much effort and some impressive successes, there has been a mixed record arising from their use 
in different applications. For example, while the Percus-Yevick (PY) equation [3,4] for a fluid of hard spheres is 
quite accurate, it proved much less successful in describing the structure of systems with longer ranged interactions 
such as the Lennard- Jones (LJ) fluid [5]. In most cases, we do not have a deep understanding of the reasons for 
a particular equation's success or failure. Part of the problem is that standard "closures" of the integral equations 
usually introduce uncontrolled approximations made mostly for mathematical convenience. Thus it is difficult to 
assess the physical consequences of the errors introduced and the kinds of interactions for which a particular equation 
is likely to be accurate. 

However, as pointed out by Stell in one of his earliest papers [6] , there is a very simple and physically suggestive 
way to interpret one of the most basic and successful of the integral equations, the PY equation for hard spheres. Stell 
noted that one can completely determine the pair correlation function h{r) of the hard sphere fluid for all distances r 
by specifying only the tail or out part of the direct correlation function c(r) (i.e., its value at separations r > d, with 
d the hard core diameter of the hard spheres). Here h and c are related by the usual Ornstein-Zernike (OZ) equation 
[1]. See Sec. (HI) below for precise definitions and further discussion. If, following OZ, one further assumes that the 
direct correlation function has the range of the potential, then its out part vanishes for hard spheres. Then the core 
or in part of c(r) for r < d can be determined directly from the OZ equation and the exact condition imposed by the 
hard core potential that h{r) = —1 for r < d. Stell showed that the resulting h{r) computed from the OZ equation is 
identical to the PY solution for hard spheres. However this simple picture directly applies only to the PY equation 
for hard spheres. 

Stell and other workers [7] generalized this idea to apply to potentials with a hard core and a longer ranged tail by 
making simple assumptions about the functional form of the out part of c(r) and solving the OZ equation subject to 
the "core condition" h{r) = —1 inside the core. The resulting mean spherical approximation (MSA) and generalized 
MSA (QMS A) equations have proved useful in a variety of applications. Madden and Rice [8] showed how these 
ideas could be apphed to systems with softer repulsive cores with their soft MSA (SMSA) equation, though the 
relationship between the original hard core condition and the treatment of soft cores, both in the initial work and in 
later derivations [1], seems (to us at least!) somewhat unclear. Most recent work on integral equation closures has 
focused attention on another function, the bridge function (see Sec. VIII C below), which is not simply related to 
the tail of c, and connections to the earlier work and the insights gained therein have often not been exploited. 
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In this paper we show how George Stell's original ideas [6] can be extended in a very natural way to describe 
more realistic systems with finite ranged soft core interactions and/or weaker and longer ranged (usually attractive) 
interactions. While some of our conclusions have been noted before, the general perspective and the formalism we 
develop is new. It gives a unified and physically suggestive way of interpreting and assessing many earlier approaches 
and ideas and suggests new and simpler approximations. The main idea is to introduce a new continuous function 
T(r) which reduces exactly to the tail of c(r) outside the (soft) core region. We show that both h{r) and c(r) depend 
only on the "out projection" of T(r): i.e., the product of the Boltzmann factor of the repulsive core potential times 
r(r). Essentially then, we have only to prescribe T outside the core, i.e., fix the tail of c, to determine h and c 
everywhere. This conclusion is rigorously true for hard cores, as noted in the original work of Stell and others [6,7]. 

We thus make direct contact with a wide class of integral equations related to the PY equation for hard spheres 
and the MSA and find in a new and more straightforward way equations related to the SMSA of Madden and Rice 
[8]. Our general approach suggests how to improve the behavior of the SMSA equation at low densities and gives 
new insights into reasons for the success of some of the most accurate integral equations, including the reference 
hypernetted chain (RHNC) equation suggested by Lado [9] and the method of Zerah and Hansen [10]. Equally 
important, many of the inherent limitations of all these methods are clarified. 



We consider here the simple case of a one component uniform fluid interacting through a spherically-symmetric 
intermolccular pair potential w{r) — uo{r) +ui(r), where uq is a harshly repulsive core potential with finite range a 
(so 7io(r) = for r > a) and ui is a longer-ranged and more slowly varying (usually attractive) potential. We will 
refer to a system with potential uq alone as the reference system and the potential ui as the perturbation potential. 
Though many of these ideas can be directly applied to fluids with long-ranged (e.g., Coulomb) forces, several new 
issues arise there that merit a more detailed discussion, and we will restrict our work here to the case where ui{r) 
goes to zero at large r faster than r^'^ . We also assume in most of the following that ui is continuous, with at least 
one continuous derivative at r = ct. Examples of a pair potential divided in this way are the separations proposed by 
Rce ct al. [11] and by Weeks et al. [12] for the LJ potential. 

The local density at a distance r away from a particle fixed at the origin in a fluid with average (number) density p 
is given by pg{r), where g{r) is the radial distribution function. In the following, we will use the notation g{r; [w]) to 
indicate the functional dependence of g{r) on the pair potential w; the subscripts will denote the reference system 
and d a hard sphere system with diameter d. Note that g{r) becomes very small in the core region r < a because 
of the repulsive core potential uq. In the special case where Uq is replaced by a hard sphere interaction Ud{r), then 
g{r; [ud + Ui]) = for all r < d. Our goal is to determine quantitatively the pair correlation function h{r) = g{r) — 1 
for the uniform fluid. Important thermodynamic and structural information are contained in h{r) and its calculation 
has been a major focus of research in the theory of liquids [1]. 



To that end most modern approaches introduce several other related functions. Probably the most fundamental 
of these is the direct correlation function c(r), defined in terms of h{r) by the Ornstein-Zernike (OZ) equation 



By iterating this equation h can be represented as a sum of chains of "direct" correlations c. For typical short 
ranged potentials, this suggests that c could be both shorter ranged than h and simpler in structure [1]. Indeed, 
Ornstein and Zernike [13] assumed that c had the range of the intermolccular potential in developing their theory of 
correlations near the critical point. While scaling theory shows that c must in fact decay as a power law r~'^ at the 
critical point, Stell and co-workers [14] have shown that very accurate results can be obtained for thermodynamic 
properties of the lattice gas surprisingly close to the critical point by assuming c is strictly the range of the potential 
and choosing its form to yield self-consistent thermodynamic predictions. Moreover, for the long ranged Coulomb 
potential, assuming that c is proportional to the potential physically incorporates the effects of screening and yields 
a nonlinear version of Debye-Hiickel theory [1] . 

We refer to the idea that c has (to a good approximation) the range of the potential as the range assumption. A 
very direct but primitive strategy for calculating h is to guess the form of the presumably simpler function c, perhaps 
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guided by the range assumption, and then determine h from the OZ equation. However, Stell's interpretation of 
the PY equation for hard spheres [6] suggests a simpler possibihty: perhaps we have to prescribe only the tail of 
c outside the range of the harshly repulsive core potential uq to determine h. We now develop a general formalism 
incorporating this idea for a system with potential w(r) = uo{r) + ui(r). 



IV. CORE AND TAIL PROJECTIONS USING CONTINUOUS FUNCTIONS 

To help us focus on the core and tail parts of functions, we note that the Boltzmann (eo) and Mayer (/o) functions 
for the harshly repulsive core potential uo{r) act very nearly as projection operators onto tail or out (r > a) and core 
OT in (r < a) subspaces respectively, since 

eo(r) = e-"^"^^^ « 0, r < a, _ ^ ^ _ ^ 1, r < a, 

^ ■' =1, r > a, •' ^ ' =0, r > a. ^ ' 

These functions exactly satisfy one property of orthogonal projectors for all r : 

-/o(r)+eo(r) = l, (3) 

and in the tail region r > a exactly satisfy the second requirement: 

- Mr) ■ eo(r) = 0. (4) 

Moreover for small r < a well inside the core, the repulsive potential uq is very large and eo essentially vanishes. 
Thus Eq. (4) also holds in this region to a very good approximation. 

However, for soft cores there is a transition region for r near a where the r.h.s. of Eq. (4) differs significantly from 
zero. Thus strictly speaking the functions — /o and cq are not true projection operators over all space. Rather they 
divide space into two parts: a tail or out part, and a core or in part. The latter is comprised of a transition region 
for r near a and an effective hard core region at smaller r. The theory for soft cores we develop works best when 
the spatial extent of the transition region is much smaller than a, as is the case for harshly repulsive interactions. 
In the special case where there is a hard core potential Ud, the width of the transition region vanishes, Eq. (4) holds 
exactly for all r, and the corresponding functions —fd and ej, are true projection operators. Our theory for soft cores 
will go over smoothly to that for hard cores in the limit of increasing steepness of the soft core potential. 

We now rewrite our correlation functions in projected form. Though our primary focus has been on the pair of 
functions h and c, both have discontinuities at r = d when there is a hard core potential Ud- It is convenient to 
introduce two new functions that remain continuous even in this limit and from which we can determine both h and 
c. One such function we will use is well known and was originally used by Stell [6] : 

t{r) = h{r) - c(r). (5) 

t is sometimes referred to as the "indirect correlation function" [15]; its continuity even when the potential has a 
hard core region is clear since it equals the convolution integral in the OZ equation (1). From this it follows that 
the first D derivatives of Hn a D-dimensional system are also continuous at r = d even for a hard core system. For 
harshly repulsive core potentials it is easy to relate c for r < a to the core part of t: to a very good approximation 
in the effective hard core region we have 

c{r)^ h{r)[l+t{r)], r <a. (6) 

This equation is exact for a hard core potential where fd and hd = —I for all r < d. 

To determine c outside the core, we now introduce a second continuous function, which we refer to as the tail 
function T{r), whose out projection eo(r)T(r) reduces exactly to the tail of c in the out region. In the core space we 
require that eo{r)T{r) correct the small errors in Eq. (6) occurring in the transition region for soft cores. Thus we 
require for all r that T(r) satisfy: 

c(r-) = /o(r-)[l + t{r)] + eo{r)T{r). (7) 
Moreover, since g = c + 1 + 1, we have, using Eqs. (3) and (7) 

gir) = eo(r)[l + t(r)] + eo(r)T(r). (8) 
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Wc have thus rewritten c and g (or h) in projected form using the new functions t and T. While special cases of 
these equations have been suggested before [6], the general utility of such a T function does not seem to have been 
realized. The most important properties of the tail function T are clear from Eqs.(7) and (8): i) it reduces exactly to 
the tail of c in the out region; ii) both h and c depend on T only through the combination eoT; iii) T is continuous 
and differentiable. 

To see that the latter holds, let us define the cavity distribution function y{r) in the usual way [1] : y{r) = 

g+/3ju(r)^^^-j Simple analysis like that mentioned above for t{r) (sec, e.g., Rcf. [4]) shows that y{r) is a well-defined 
continuous function of r with several continuous derivatives even when w itself has a hard core region or other 
discontinuities. Using Eq. (8) we immediately get that 

y{r) = [l + t{r)+T{r)]/e,{r). (9) 

Here ei(r) = e~^"^'^'^K Since y{r) and t{r) are continuous and differentiable and the perturbation tail function 
ei(r) can be constructed to be continuous and differentiable even across a hard core region, it follows that T(r) is 
continuous and differentiable ^. When the potential has a hard core, Eq. (9) can alternatively be used to define T(r) 
for all r in terms of the more familiar functions y, t, and ei. 

V. BASIC RESULT 

Now wc can refine the primitive strategy of guessing c and using the OZ equation to calculate h, by reexpressing 
everything in terms of t and T. See the Appendix for numerical details. In principle, if we prescribe T{r) for all 
r then t{r) can be completely determined from the modified OZ equation. However, we see from Eqs. (7) and (8) 
that since both g and c (and hence also t) depend only on cqT, the results are very insensitive to any errors we 
make in prescribing T in the core space r < a. This is obvious in the effective hard core region where cq essentially 
vanishes. In the narrow transition region, since T is continuous and differentiable, its values there can be accurately 
determined by extrapolation from those for r ^ a. In effect then we only have to prescribe the out part of T, i.e., the 
tail of c, to determine both h and c everywhere. This generalizes Stell's argument [6] for the hard core PY equation. 
In the Appendix we introduce a new and very efficient variational method that allows us to determine numerically 
both h and c from the OZ equation given some approximation for the out part of T{r). This will allow us to find 
accurate solutions to many standard integral equations in a very simple way. 

Note from Eq. (9) that the tail of c is not sufficient to determine y{r). Its values for small r in the effective hard 
core region depend directly on T(r) there and we cannot expect that extrapolation from the out part of T alone will 
give accurate results for T{r) well inside the core. From this perspective, the calculation of y(r) (and other closely 
related functions such as the bridge function B{r) = lnt/(r) — t{r) discussed below in Sec. VIII C) is a much more 
difficult problem, requiring the accurate determination of both the out and core parts of T(r). Fortunately the latter 
problem does not have to be solved to find accurate results for h and c. This point was emphasized by Stell for the 
hard sphere system [6] , and we see it holds true much more generally. 

VI. RELATION TO PREVIOUS WORK 

Stell's original work [6] was designed to provide information about the PY equation for a system with the general 
pair potential w{r). To that end, he introduced a set of equations very similar in form to Eqs. (7), (8), and (9), but 
with the crucial difference that the Boltzmann and Mayer functions e and / for the full potential w appear, where 

e(r) = e-f"^^'^ = eo(r)ei(r) ; /(r) = e-'^-M - 1 = /o(r) + eo(r)/i(r) . (10) 

Here /i(r) = e~^^^^^^ — 1. Note that / has the range of the full potential and — / and e no longer approximate 
projection operators onto core and tail regions. Stell's equations can be written as 



^More generally, we can exploit the fact that y and t have at least 2 continuous derivatives for D = 3, to relate the behavior 
of low order derivatives of T to those of ei . This could be used to give a more accurate extrapolation of T into the transition 
region. 
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c(r) = /(r)[l + t(r)]+e(r)rf(r), 



(11) 



g{r) = e(r)[l + t{r)] + e{r)d{r) 



(12) 



y{r) = l+t{r)+d{r). 



(13) 



Eq. (13) can be taken as the definition of the function d{r) (we use Stell's notation; this should not be confused with 
the hard sphere diameter). Despite the superficial similarity of these equations to our Eqs. (7), (8), and (9), d{r) in 
general has very different properties than our analogous fmiction T{r). In particular, d{r) docs not reduce to the tail 
of c in the out region and is likely to have a more complicated oscillatory structure. The main utility of Eqs. (11), 
(12), and (13) is in analyzing the PY equation: Stell was able to show that the usual formulation of the PY equation 
for a general potential results from the approximation d{r) = 0. Unfortunately there is little reason to believe this 
approximation is generally accurate. 

However, in the special case of hard core interactions where w{r) = Ud{r) , Eqs. (11), (12), and (13) reduce to our 
Eqs. (7), (8), and (9), and dd{r) — Td{r). The approximation dd{r) = in the out region for hard spheres then can be 
motivated by an application of the range assumption for the tail of c. This assumption alone is enough to determine 
the accurate PY solution for hd(r). The range ansatz dd(r) = for r > d is exact in one dimension {D = 1) and 
hence yields the exact hd{r). In D = 3, the first errors in (r) show up at 0{p'^) in a density expansion. Overall 
h^^ (r) remains remarkably close to the results of computer simulations even at higher densities, with small errors 
most noticeable near contact and at the first minimum for densities near the fluid-solid transition [1]. As noted by 
Stell [6], all that is required to calculate hd{r) in general is an expression for dd{r) in the out region. Essentially 
exact results for hd{r) can be obtained from the generalized MSA (GMSA) of Waisman and Lebowitz [16], which 
assumes the existence of a small short-ranged (Yukawa-like) tail in Cd{r) for r > d. Parameters in the tail are chosen 
so that hd gives results for the pressure and compressibility that fit simulation data. The basic picture suggested by 
the range assumption that the tail of Cd has a simple structure and is small and much shorter ranged than ha seems 
to be well established. 

Stell [6] also noted that the extrapolation of the PY approximation dd{r) = deep into the core space is a separate 
and much less accurate approximation. For example, the resulting PY expression for yd{r) given by Eq. (13) with 
ddir) = for all r < d can have large errors at small r for D = 1 even though the PY result for hd{r) is exact. (While 
dd{r) is continuous and difi^erentiable at r = d higher derivatives are discontinuous, leading to a large positive value 
at small r for the exact dd{r) at high density.) This strongly suggests that the calculation of hd and yd should be 
logically separated [17]. Of course, ya is an interesting function and additional properties like the chemical potential 
can be obtained from it [1]. However, a focus on hd and Cd alone permits a very simple theory, and one can use 
results for the pressure and compressibility from gd{r) and Cd{r) and thermodynamic relations to calculate other 
thermodynamic properties. In particular, in this approac;h the chemical potential should be calculated by integrating 
the pressure, and not from the very inaccurate value for t/d(0) given by extrapolating dd{r) = deep into the core 
space. By introducing the tail function T(r) and the system of equations (7), (8), and (9), we have been able to 
extend these important ideas of Stell for hard sphere systems [6] to systems with more general interactions. 



We now describe some general properties of T{r). Using Eqs. (7), (8), and (9), this can be rewritten exactly as 



explicitly showing that T reduces to the tail of the direct correlation function in the out region, but has a different 
form in the core region. To focus on the changes induced by the perturbation potential ui, it is useful to define the 
excess quantities: 



where Tq is the exact T function for the reference system, with similar deflnitions for other excess functions such 
as Ah and Ac. According to the range ansatz Tq is zero in the out region, and we expect that the exact Tg will in 
general be small and vanish rapidly at larger r outside the core. Thus in the out region T(r) « AT(r), and is mainly 
determined by the potential tail ui{r). 



VII. GENERAL PROPERTIES OF THE TAIL FUNCTION 



r(r) = c(r) - /o(r)ei(r)y(r), 



(14) 



AT(r) = T(r) - To(r), 



(15) 
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Based on an analysis by Stell [18], it is generally believed that away from the critical point the asymptotic form 
of c(r) at large r is 



c(r) ~ -/3ui(r). (16) 

For system with a weak and slowly varying potential tail vi that goes smoothly to zero at large r this is consistent 
with the idea that the OZ equation should reduce to linear response theory far from the core region. Here /3 is the 
inverse of Boltzmann's constant times the temperature. Thus we expect AT(r) —f3ui (r) far from the core. 

At very low density p graphical expansion methods show that the exact form of c(r) for interaction potentials 
going to zero faster than can be written as: 

c{r)=f{r)[l+pA{r)]+0{p% (17) 

where 



A(ri2) = J dr3/(ri3)/(r32). (18) 

Note that the range assumption for c is rigorously true at low density. Similarly it is easy to show that 

t{r)=pA{r)+Oip% (19) 

y(r) = l + pA(r)+0(p2), (20) 

and 

T{r)=,h{r)[l + pA{r)]+0{p'). (21) 
It follows from Eq. (21) that To(r) = + 0{p^). 



VIII. CLOSURES AND THE TAIL FUNCTION 



Most integral equation theories for h{r) are based on the idea of a closure [1]: a second relation between h and c 
which, when combined with the OZ equation, allows one to solve for the values of h and c. However most closures 
are expressed in terms of more complicated functions like y{r) or B{r) and their form is usually determined by 
mathematical considerations. See, e.g., Sec. (VIHE) below. The above results show that to calculate h{r) we can 
focus on the simpler projected function eo(r)T(r), determined essentially only by the tail of c(r). An exact choice 
will yield an exact h and approximate choices can be motivated by the range ansatz and the general supposition that 
the tail of c has a simple structiire. As discussed in the Appendix, we can also exploit the relatively simple nature 
of the out part of T{r) in the numerical solution of the resulting integral equations. Other standard closures can be 
reinterpreted and sometimes simplified by looking at their predictions for the tail of c. 



A. Soft Mean Spherical Approximation 

Probably the simplest such prediction directly yields the SMS A integral equation [8]. The SMSA assumes that 
the limiting linear response value for the tail of c given in Eq. (16) holds for all r in the out region. Thus we set 

eo(r)r^^^-4(^) = eo(r)[-/3ui(r)] (22) 

in Eqs. (7) and (8). In the out region we have Tq^^^ = and /\j;'SMSA _ 'j'j-^g resulting expressions for 

h and c can easily be shown to be equivalent to the original SMSA results, which were written in a different form. 
If ui = then the SMSA reduces to the PY equation for the reference system. The approximation T^^SA = g in 
the out region again can be motivated by the range assumption. When uq is replaced by a hard core potential Ud 
then Eqs. (7) and (8) with Eq. (22) reduce to the original hard core MSA. This derivation and interpretation of the 
SMSA and its relation to the MSA seems much simpler than that found in previous work. 
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One way to improve the SMSA is to improve its description of repulsive forces. Equation (22) sets Ti^'^^^A — g in 
the out region. If a more accurate expression for Tq is known this could be used along with the MSA approximation 
^rpSMSA _ in the r.h.s. of Eq. (22). For hard cores the GMSA [16] should give a very accurate expression 

for Td{r). Its use in the r.h.s. of Eq. (22) for a system with potential w = Ud + ui would yield a theory essentially 
equivalent to the optimized random phase (ORPA) theory of Andersen and Chandler [19], where exact hard sphere 
correlation functions are supposed to be used along with a MSA treatment of wi . 

The SMSA gives rather accurate results for the high density LJ fluid and correctly describes the qualitative changes 
in Ah = h — ho induced by wi. However, it is much less accurate at low densities. This can be understood since 
Eq. (22) does not reduce to the exact result, Eq. (21), at low densities. An improved theory would result from 
approximations for T(r) that interpolate between the exact low density limit, Eq. (21), and Eq. (22) at high density. 
We will describe several such theories below. 



B. PY and HNC Equations 

Other integral equation closures can be reexpressed in terms of their predictions for the out part of T. In many 
cases this can give us insights into their strengths and weaknesses. For example, by rewriting the standard expression 
gHNC _ exp(— /3i« + 1) given by the hypernetted chain (HNC) equation [1] in the projected form of Eq. (8), we find 
that the HNC closure predicts 

T"^^ = exp(-/3ui +t)-{l+ 1). (23) 

This agrees with the exact Eq. (21) at low density. However, when applied to the reference system, Eq. (23) predicts 
that T^^'-^ — exp(to) ^ (1 + io)- Since is large and oscillatory at higher density in the out region, this strongly 
violates the range assumption. Indeed the HNC equation gives very poor results for a dense hard sphere system. 
Experience has shown that the HNC closure does a much better job of describing slowly varying interactions, and 
for systems with long-ranged Coulomb forces it is often the theory of choice [1]. As discussed below one of the most 
accurate integral equation theories, the RHNC theory [9], combines a HNC treatment of the more slowly varying 
potential ui along with an (in principle) exact treatment of reference system correlations. 

The PY closure for the reference system incorporates the range assumption and gives a much better description 
of reference system correlations than does the HNC. However, for the full system it predicts for the out part of T: 

T'''' = h{\ + t). (24) 

This again agrees with Eq. (21) at low density. However at higher density the oscillations in t and the strong nonlinear 
dependence on the perturbation potential will yield a larger and more oscillatory tail for c than suggested by the 
SMSA in Eq. (22). In practice the simple linear response form of the SMSA gives much more accurate results at 
high density [8]. 



C. Bridge Function 

Most recent integral equation closures focus attention on another continuous and differentiable function, the bridge 
function B(r), which can be defined formally as [1] 

B{r)=\ny{r) -t{r). (25) 

Thus g{r) = cxp[— /3w(r) + t{r) + B(r)\. B(r) represents the sum of a well-defined set of Mayer cluster diagrams, 
and the HNC equation results from the approximation B{r) = 0. B plays a role analogous to our function T in 
generating closures, and we shall see that some of its relevant properties can be understood more easily from those 
of T. Thus one can represent h, c, and y in terms of the pair of functions B and t. If B is specified by some closure 
ansatz, then these functions can be calculated using the OZ equation. 

Alternatively, using Eqs. (7), (8), and (9), we can exactly express B in terms oft and T: 

B{r) = ln[l + t{r) + T{r)] - [t{r) - I3ui{r)]. (26) 

Thus B depends on T itself rather than the projected function cqT, and in that sense is a more complicated function 
than h or c. Indeed determining its form, particularly inside the core, has proved a very difficult challenge both for 



7 



theory and simulation, and definitive results are still not known [20]. However, since the out part of B in Eq. (26) 
can determine the out part of T, we can effectively concentrate only on the out part of B if we restrict ourselves to 
theories for h and c. 

In general, the out part of B has a rather complicated oscillatory structure. For example, for the reference system 
we have exactly in the out region, using the definition of t, and the equality of the tails of c and T, 

Bo{r) =\n[l + ho{r)]-ho{r)+To{r), r > a. (27) 

Since the exact To is almost certainly small and very short ranged, as suggested by the range ansatz and the success of 
the PY equation for repulsive forces, Bq will have longer ranged oscillations determined by those of the pair correlation 
function ho- Setting Tq = in Eq. (27) yields the PY expression for the reference system bridge diagrams. 

However, in many cases the oscillatory tail of B{r) for the full system seems to depend only weakly on the 
perturbation potential ui{r), so that B{r) « Bo{r). This idea has been called the universality of the bridge function 
[21], with Bo often approximated by B^, the bridge function of an appropriately chosen hard sphere system ^. The 
following argument gives some insight into why this could be a reasonable approximation for the out part of B. 
Analogous to Eq. (27), we have exactly 

B{r)=\n[l + h{r)]-h{r) + [T{r)+(3ui{r)], r > a. (28) 

At high density, the structure is dominated by repulsive forces for systems with short-ranged interactions [12] and it 
is a fairly good approximation to set h{r) « ho{r) ("universality" of the correlation functions!) Moreover the success 
of the SMSA suggests that T(r) « — /3ui(r) and To(r) w are also reasonable approximations in the out region. 
Then Eqs. (27) and (28) yield B{r) « Bo{r) in the out region. Note that this result is exact at low density since 
B = + 0{p^). Thus for this class of systems, we can arrive at the idea of approximate bridge function universality 
outside the core using the more physically transparent arguments of the SMSA. Differences in the results for the 
two theories should be small at high density. It can be seen using the general expression for B in Eq. (26) that 
these arguments do not hold for the core part of B and we see no reason to expect any such "universality" at higher 
densities there. 



D. RHNC Equation 



Alternatively, if we assume it is a good approximation to set B{r) « -Bo(r) in the out region, then for systems 
where h{r) « ho{r) we have T(r) w —0ui{r) from Eqs. (27) and (28), which is the SMSA closure. At low density 
h(r) « ho{r) is not accurate, and the true T(r) must differ significantly from the SMSA prediction. Indeed using 
the exact low density forms for h and ho along with B{r) = Bo{r) in Eqs. (27) and (28) yields the exact low density 
form for T given in Eq. (21). Thus a theory incorporating B{r) w Bo{r) in the out region will give exact results for 
h at low density and should give results at high density close to those of the accurate SMSA. 

This is what is done in the RHNC theory of Lado [9] , and overall this is one of the most successful integral equation 
methods known. The standard RHNC closure can be written as 

^flHJVC(^) = exp[-/3w(r) + t{r) + Bo{r)], (29) 

thus replacing the exact bridge function Bhy Bq. To describe its predictions in terms of T, it is convenient to consider 
excess functions like that defined in Eq. (15). We find 

eo(r)AT«^^^(r) = go{r){exp[-(3ui{r) + At{r)] - 1} - eo{r)At{r). (30) 



^In applications to tfie RHNC equation, discussed in Sec. VIII D, the hard sphere diameter d is often taken as a parameter 
that can be varied to achieve more consistent thermodynamic predictions from the full system's correlation functions. However, 
for the systems we consider here with short-ranged interactions, it seems more realistic to fit d to properties of the reference 
system using, say, the blip function expansion [1]. For more accuracy, one can directly approximate Bo using various closures 
that accurately describe soft repulsive systems, as suggested in Ref. [11]. Coulomb systems with strong long-ranged repulsive 
and attractive forces require special treatment, and can have correlation functions differing considerably from those of the 
reference system. In such cases, choosing d to represent some effective hard core diameter for the full system may be a 
reasonable first approximation. 
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Note that we only require accurate values for .90 (r) and not for Bo{r) well inside the core to determine this fundamental 
quantity in our approach ^ . A numerical solution can be found using the general variational method described in the 
Appendix. 

To examine the relation between the RHNC and the SMSA more quantitatively, let us define 

AT(r) = -/3«i(r)+C(r). (31) 
For the SMSA ^smsa^^-^ = in the out region. We find that in the out region ^^^hnc^j.^ ^g^^ written exactly as 

$^™(r) = A/i(r)-ln[A/i(r)/5o(r-) + l], r > a. (32) 

This agrees with exact results from Eqs. (21) at low density and corrects the poor behavior of the SMSA there. At 
higher density, ^^^^'^ represents an additional oscillatory component in the tail of T when compared to the SMSA.. 
However, when Ah is small, as is generally the case at high density for the systems we consider, then ^^hnc jg 
small (with ^^^^c- vanishing whenever A/i(r) = 0). Thus AT^^^*^ » ^j^smsa ^ -(3^ in the out region at high 
density, as argued above. 

E. Unique Function Ansatz 

Several workers have tried to find more accurate expressions for B{r) by assuming it is some unique local function 
[22] of t{r), as suggested by several approximate closures that gave good results for systems with short ranged 
repulsive interactions [23] . LLano-Restrepo and Chapman (LC) showed for systems with an attractive potential tail 
ui that this assumption was generally inaccurate at small r in the core region and also was inaccurate at high density 
outside the core [24] . They proposed that there could exist some "renormalized" function t{r) involving m such that 
B is a local function of t. They found that the choice 

t{r)=t{r)-l3ur{r) (33) 

gave accurate results at high density in the out region for the LJ fluid. This is precisely what would have been 
suggested by applying the SMSA closure T{r) = —0ui{r) to the exact Eq. (26) in the out region. 

However, the SMSA approximation for T is not accurate well inside the core space and indeed the renormalized 
function gave poor results there. Moreover the SMSA approximation for T is inaccurate at low density where the 
exact T reduces to /i . Indeed this shows that the local function ansatz for B cannot in general be correct even outside 
the core. Duh and Haymet [20] and Duh and Henderson [25] have proposed different density dependent separations 
of the total potential: w{r) = Uo{r;p) + ui{r;p), with "reference" (uq) and "perturbation" (ui) parts chosen such 
that Eq.(33), now defined with ui, could give more accurate results for _B as a local function of t, even well inside the 
core where LC's original suggestion most noticeably failed. It is clear from Eq. (26) that the unique function ansatz 
can give exact results at low density only if the perturbation Ui{r; p) vanishes as p — > since then T ^ Tq = 0, as 
shown in Refs. [20] and [25]. Assessing the nature of errors induced by the unique function approximation in general 
remains a very difficult task. For our purposes here it seems simpler and more direct to retain the original physically 
motivated separation and focus instead on the out part of T, whose density dependence is such that T reduces to /i 
at low density while approximating — /3ui at high density. 

IX. CLOSURES SATISFYING CONSISTENCY CONDITIONS 

A natural idea is to consider more general density-dependent expressions for T that can vary between these limits, as 

suggested by the RHNC equation. Parameters in the interpolation function can be chosen to fit simulation data or to 
satisfy various thermodynamic consistency conditions (Maxwell relations and sum rules) which the exact correlation 
functions must obey. We first discuss one of the most successful integral equation approaches, the method of Zerah 
and Hansen (ZH) from this perspective [10], and then introduce a new and simplified method which implements this 
idea in a very direct fashion. Results seem very promising. Contact is also made with very recent work by Stell and 
coworkers [26]. 



^An even simpler approximation in the spirit of the RHNC equation suggests itself, where go{r) in Eqs. (30) and (32) is 
replaced by eo(r). We expect this to have essentially the same behavior at both high and low density. 
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A. HMSA Equation. 



ZH introduced a generalized "HMSA" closure that interpolates nonlinearly between the SMSA closure at small r 
and the HNC closure at large r, with a parameter in the interpolation function chosen to give consistent results for the 
pressure computed from the virial and compressibility formulas [10]. The choice of the HNC theory at large distances 
was motivated by its superior behavior for systems with long-ranged forces. The ZH closure can be rewritten as the 
following expression for T in the out region: 

rj^zH.. _ exp{f„(r)[t(r) - (3u,{r)]} - 1 - F^{r)t{r) 

where Fa{r) is an r-dependent interpolation function, 

F„(r) = l-exp(-r/ra), (35) 

and Va a fitting parameter chosen to achieve thermodynamic consistency. For Fq — !■ ( i.e. for r/r^ 0) Eq. (34) 
reduces to the SMSA closure —(3ui and for Fa ^ 1 { i.e. for r/r^ oo) Eq. (34) reduces to the HNC closure, 
Eq. (23), though Eq. (35) implies a rather slow transition between these limits for physically relevant values of r^. 

For the systems we consider here with short-ranged interactions, the important feature of Eq. (34) is not the 
behavior of the HNC equation at large distances but the fact that at low densities T^^'-' reduces to the exact result, 
Eq. (21). ZH found numerically for the LJ fluid that decreased as the density tended to zero, so the HNC closure 
is effectively used at all relevant r at very low density. At higher density increases, thus mixing in more and more 
of the SMSA expression. For example, near the triple point (at a reduced density of 0.85 and a reduced temperature 
of 0.786) ZH found that = 6.25a- [10]. The ZH interpolation scheme provides a mechanism by which one can go 
between these limits as the density changes while maintaining enough flexibility in the shape of T outside the core 
that thermodynamic consistency for the pressure can be achieved. 



B. Tail Interpolation Method 

Both the ZH equation and the RHNC equation discussed above in Sec.VHID give accurate results at both high 
and low densities by considering some rather complicated density dependent expressions for the out part of T, which 
in particular involve t{r). See Eqs. (34) and (30). The variational method discussed in the Appendix can be used to 
solve the OZ equation when the out part of T is a known function of r, as is the case for the SMSA approximation. 
Because of the appearance of the initially unknown function t{r) in the ZH and RHNC expressions for T(r), we 
cannot use this method alone to solve these equations. However, by making an initial guess for the out part of T, 
we can iterate until the value of the out part of T does not change. This method combines the standard Picard 
iteration scheme for the hopefully slowly varying out part of T with the efficient variational method for the core parts 
of functions. While we have found that this method generally works quite well, it still requires much more computer 
time than does the direct variational solution of the OZ equation with a known out part of T. Moreover because of the 
complicated nonlinear nature of the self-consistency condition and the direct interplay between possible oscillations 
in t and T in the out region it is not clear that self-consistent solutions can always be found for physically relevant 
states. Indeed the RHNC equation fails in a quite peculiar way [27] close to the vapor-liquid coexistence region. 

We now introduce a new method, which we call the tail interpolation method (TIM), that implements the idea of 
a density dependent interpolation involving /i and — very directly, while using a very simple {t independent) 
expression for the out part of T. We assume the out part of AT can be written as: 

^TTIM ) ^ a/i (r ) + ( 1 - a) [-/3wi (r )] , (36) 

where a is a (temperature and density dependent) parameter that is chosen so that consistent results for two different 
routes to one particular thermodynamic quantity are obtained. (To obtain the full T the out part of Tq should be 
added to Eq. (36); often the SMSA-PY approximation Tg = gives sufficient accuracy.) Note that the presumably 
exact asymptotic form for the tail of c(r) given in Eq. (16) is maintained for any choice of a, and a is not required 
to lie between zero and one. In general, varying a allows us to change the shape of the tail of c at intermediate 
distances while maintaining the proper asymptotic form, and we use this freedom to achieve partial thermodynamic 
self-consistency. At low density a = 1 and, given the relative accuracy of the SMSA, we expect that at high densities 
smaller values of a will be found. 
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In this paper wc impose consistency between the virial and compressibihty routes to the isothermal compressibility. 
Belloni [28] has shown that this can be implemented very efficiently by differentiating the OZ equation, and our 
variational method can be easily extended to this case. We have not yet examined in any detail the merits of 
this choice over other possible consistency conditions. Indeed the SMSA usually gives rather poor results both for 
the virial pressure and for the compressibility [30], and the energy route is typically used to give more accurate 
thermodynamic results [26]. It is easy to derive a variational method to impose thermodynamic consistency from 
the energy route and we suspect this will give even better results. However, in this initial study we have imposed 
consistency on the isothermal compressibility to see whether self-consistency using the very simple expression for 
T(r) given in Eq. (36) can improve on the rather poor performance of the SMSA for this quantity. The preliminary 
data we report in the next section illustrates the basic concept and suggests that further work is indeed merited. 

X. NUMERICAL RESULTS 

We test our approach on two well-studied systems: the hard core Yukawa fluid (HCYF) and the LJ fluid. The 
HCYF has been the focus of recent theoretical work [26] and represents a system where errors from the treatment of 
soft cores do not arise, while the LJ fluid is a typical soft core system. 

A. HCYF 

The interaction potential in the HCYF is given by: 

«;^c^^^(r-)=u,(r)+e^^^, (37) 

where d is the hard sphere diameter. Wc choose z = 1.8/d, which corresponds to a well-studied system [29,30]. 
We have solved the TIM equations using the variational method described below in the Appendix. For greater 
accuracy in treating the hard sphere correlations at high density, we have included a GMSA like expression for 
in the out region, as described in the Appendix. Only preliminary results are reported here. In Fig. 1 we give 
values for the compressibility factor pP/p compared to the results of a MD simulation study [30]. We emphasize 
that the compressibility factor has been calculated directly from the virial formula for pressure and not obtained 
through thermodynamic relations from the energy route, as is usually done in ORPA and MSA approaches for greater 
accuracy. In the inset to Fig. 1 we present the dependence of the TIM self-consistency parameter a on temperature 
and density. Isotherms T* = 2.0 and T* = 1.5 are supercritical, and T* = 1.0 is subcritical ^. At low densities a 
approaches the exact limit a = 1, while at higher densities a becomes smaller though differing from zero (the MSA 
limit). The behavior at intermediate densities where a reaches a maximum is interesting and was not anticipated by 
us. The behavior of the TIM very near the critical point and spinodal lines has not been examined. 

To test the accuracy of the correlation functions predicted by the TIM, we compare them to the results of new 
MC simulations we have carried out ^. In Fig. 2 we show h{r) given by the TIM, the ORPA, and an even simpler 
self-consistent approach (SC2) very similar to that used by Stell and coworkers [26], where T'^'^^(r) = a[— /3wi(r)], 
with a is chosen to satisfy self-consistency of the virial and compressibility routes to the compressibility. Since the 
SC2 tail does not have enough flexibility to reduce to /i at low density, we expect that its correlation functions may 
be less accurate there. The results show the relatively inaccuracy of the ORPA at intermediate and low densities, 
with best results seen at high density. The SC2 approach, while giving accurate self-consistent thermodynamics, 
yields less accurate correlation functions at low densities, as expected. 



''We use reduced units: p* = pd^, T* = ksT/e., etc. 

^This is a standard NVT-ensemble MC simulation, with the number of particles ranging from 128 to 432, depending on 
density. The pair potential has been cut and shifted at rcut/d = 3. 
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B. LJ fluid 



We have also solved the TIM equations for the LJ fluid for a few states, using the WCA separation of the pair 
potential [12]. For the relatively low density states we study here, the SMSA-PY approximation for the reference 
system Tq = gives sufficient accuracy. In Fig. 3 we compare predictions of TIM and SMSA approximations to MD 
simulations results [31]. The states shown correspond to low and moderate densities at about the critical temperature. 
Again the TIM approach gives notable improvement over the SMSA theory, especially at low densities. 

XI. FINAL REMARKS 

Many issues in the theory of integral equations for fluid structure can be profltably analyzed and interpreted in 
terms of predictions for the out part of the tail function T(r), i.e., the tail of the direct correlation function. In this 
paper we have only considered a single component uniform fluid with short ranged interactions. Here the simplest 
possible MSA linear response approximation relating the tail of T(r) to the perturbation potential ui{r) immediately 
yields the SMSA theory. For systems with harshly repulsive forces only, the SMSA reduces to the successful PY 
theory. For systems with a weak potential tail wi(r), the SMSA gives rather accurate results at high density but fails 
at low density. The behavior at low densities can be greatly improved by introducing a density dependence into the 
tail of T{r) such that it reduces to the exact low density result of Eq. (21), as is effectively done in the RHNC and 
HMSA equations. We introduced here a new self-consistent (TIM) method that incorporates this idea in a much 
simpler form, and the preliminary results for the LJ fluid and the HCYF appear promising. 

The success of all these methods at high density arises from the fact that for the systems considered attractive 
forces have only a relatively small effect on the liquid structure, so that h(r) ~ ho{r) is a fairly good approximation. 
Put another way, the density fluctuations can be well described by a simple Gaussian theory [37] . When this is not 
true, as is the case for nonuniform liquids [39], particularly in cases of wetting and drying phenomena, the natural 
(singlet) generalizations of all these integral equation methods fail [38] . We simply do not know a good enough guess 
for the tail of T(r) in cases where attractive forces induce such signiflcant structural changes. New approaches based 
on a self-consistent mean field treatment of the attractive interactions appear more promising here [39] . 

A more severe test of these ideas and of the utility of integral equations in general for uniform fluids is in applications 
to systems with long-ranged Coulomb interactions. Here one must deal with fluid mixtures with strong and long- 
ranged attractive and repulsive interactions, and the correlation functions differ in significant ways from those of 
any reference system with short-ranged forces. Nevertheless, characteristic properties of systems with long-ranged 
forces, such as the Stillinger-Lovett moment conditions are very naturally expressed in terms of the tail of the direct 
correlation function [1]. The RHNC and HMSA approximations have often proved useful here, and even the simple 
SMSA captures Debye screening, perhaps the most fundamental feature of the long-ranged force problem. We hope 
that some the ideas presented here for the T function can be extended to long-ranged force systems to provide a 
more intuitive understanding of the strengths and weaknesses of existing integral equation approaches, and aid in 
the development of new and simpler approximations. 
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APPENDIX A: VARIATIONAL METHOD 
1. Fixed Tail Function 

We can look on the OZ equation (1) as indirectly relating the continuous functions T(r) and t{r). Thus, given 
T(r) we can in principle solve for t{r) and then determine h{r) and c(r) from Eqs. (7) and (8). We first consider the 
simplest case, exemplified by the SMSA, where the out part of T(r) is a fixed prescribed function, independent of 
other correlation functions and/or the density. Then we show how to generalize this approach for an arbitrary choice 
of T(r), which can be coupled to other correlation functions such as t{r). In the latter application our approach 
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represents a new way to solve standard integral equations, and we believe it offers some notable advantages over 
conventional methods. 

It is easy to rewrite the OZ equation (1) in terms of t and c. Taking Fourier transforms we have 

where c{k) denotes the Fourier transform of c(r). As noted above, only the "out projection" eo(r)T(r) is actually 
relevant for h and c. Given this, Eq. (7) shows that wc need to fix only the "core projection" /o(r)f(r) to determine 
c(r) for all r. In principle, t{r) can then be determined everywhere from the modified OZ equation (Al). A proper 
self-consistent choice for t{r) inside the core must yield the same functional form when it is computed indirectly 
using the OZ equation (Al). This requirement can be formulated very efficiently in terms of a variational procedure. 

In the following analysis eo(r)T(r) is held constant and variations in all functions are generated solely by variations 
in t{r) restricted to the core region r < a. According to Eq. (7), variations of c(r) and t{r) then are linearly related: 

6c{r) ^ Mr)6t{r). (A2) 

To arrive at the proper variational functional, we first formally integrate the r.h.s. of Eq. (Al) with respect to c, thus 
arriving at a functional 

'^oz = J {p'c\k)/2 + pc{k) + ln[l - pc{k)]}d\<i , (A3) 

whose general variation with respect to c can be simplified using the modified OZ equation (Al): 



(27r)3 / 



P c{k) + p ■ 



Sc{k)dk 



- pc{k)_ 

i{k)dc{k)dk . (A4) 



(27r)3 

Using Parseval's formula, and considering the special variation of c given by Eq. (A2), we have 

S^oz = P^ J fo{r)t{r)St{r)dr , (A5) 

which expresses the result in terms of the imposed variation of t inside the core. In Eq. (A5), t{r) satisfies the OZ 
equation (Al). We now consider a second functional of t : 



^dir = -y / h{r)t\r)dv (A6) 



whose variation directly gives the negative of the r.h.s. of Eq. (A5). Thus by construction, the functional 

$ EE $oZ + ^d^r (A7) 

obtained by adding Eqs. (A3) and (A6) is stationary (and reaches its minimum) when the proper self-consistent 
value for t(r) inside the core is used. 

To implement this variational method, we expand the core part of t{r) iov r < a in terms of Legendre polynomials, 
orthogonal on [0,0-]: 

n 

t(r) =^a„F„(r), r<a. (A8) 

i=l 

We choose values of the coefficients a„ to minimize the functional $ in Eq. (A7). We have used Powell's quadratically 
convergent method to implement the minimization procedure [32]. If needed, one could improve this step of the 
calculation by using conjugate gradient methods. In practice, our implementation is very efficient, and t is smooth 
enough that it is generally sufficient to use n « 5 to get highly accurate results. For example, for hard core systems 
— {l+t) = c inside the core, and the exact solution of PY equation for hard spheres gives a c that is a cubic polynomial 
[1]. Fast Fourier Transform methods [32] are used in evaluating Eq. (A3). 
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One important general feature of our method is that we expand the smooth function t inside the core rather 
than c. As discussed above, for hard core systems these two procedures are equivalent, and our method then reduces 
exactly to the variational method Andersen and Chandler [19] used to solve the hard core MSA and ORPA equations. 
However, for soft core systems, while c(r) is simply related to t{r) well inside the core, it changes rapidly in the 
narrow transition region, close to r w ct. Higher order polynomials are required to describe this rapid localized 
variation of c accurately. This problem becomes more and more severe with increasing steepness of the reference 
system potential. However t changes smoothly and slowly even in the transition region. This is illustrated in Fig. 4, 
where we show t{r) and c(r) both for the full LJ system, calculated using the self-consistent TIM method discussed 
in Sec. (IX B), and for the LJ reference system using the SMS A (PY) approximation. The same qualitative features 
are seen both at lower and higher densities and using other accurate closures. Previous variational methods proposed 
for soft core systems [33,34] have either expanded c(r) inside the core or h{r) itself. 

2. Arbitrary Tail Function 

To solve integral equations for an arbitrary prescription for T(r) containing initially unknown functions like t{r) 
(see, e.g., Eqs. (24), (23), and (30) for the PY, HNC and RHNC equations) we can combine the variational technique 
with an iterative method. First, we make an initial guess T^^^r) for the out part of T, and solve the variational 
problem as above for this fixed choice. This will yield new values for t and other correlation functions. Then, the 
next approximation T'^^)(r) can be determined for a given closure and the obtained correlation functions. The new 
approximation is used in the next variational step and this iteration procedure is repeated until convergence of the 
tail of T(r) is obtained. Again one could replace the iterative steps by more sophisticated methods [35]. However, 
since the out part of T has a relatively simple structure, we have found the simple iterative method works quite well 
in all cases we have tested. 

3. Thermodynamic self-consistency 

By introducing a dependence of T(r) on one free parameter we can ensure (partial) self-consistency of thermody- 
namic properties. Here we extend our variational method to impose consistency between the virial and compressibility 
routes to the isothermal compressibility: 

X?(p,/3)=X?(p,/3). (A9) 

Here 



Xt(p,/3) 



To evaluate Xt(P' I^) need an efficient way to calculate dg{r)/dp. Just as we did for g{r) we can use a variational 
method. 

To simplify notation, let us denote the density derivative of a function f{r;p) as fp{r) = df{r;p)/dp. The OZ 
equation (Al) and relation (7) can be directly differentiated: 



\k)+ pc{k)[2-pc {k)]c'p{k) 
[1 - pc(fc)]2 



4(r-)=/o(r)t;(r)+eo(r)T;(r), (A13) 

5;(r)=t;(r)-c;(r). (A14) 

To solve (A12) and (A13) one needs to know the function Tp(r). In the simplest approach, this derivative is 
neglected, because its density dependence is usually very weak. This has been referred to as local consistency [28] . 
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Alternatively, one can first calculate T(r) with Tp{r) set to zero and then compute its density derivative by a finite 
difference method by evaluating T{r) at slightly different densities. The results can be plugged back into (A13) and 
the equations iterated until convergence to globally consistent results, like those in Ref. [26] , are found. 

In calculations of correlation functions of the HCYF we have introduced a small correction to the out part of T(r), 
which corresponds to a proper description of the pure hard-sphere system (the limiting case of HCYF as /3 ^ 0). 
This correction has the usual GMSA-like Yukawa form [16] 

Td{r) = Kd[p) exp[-2d(p)(r - d)]/r, (A15) 

with Kd{p) and Zd{p) chosen to satisfy consistency of the virial and compressibility pressure and the Carnahan- 
Starling equation of state. In practice, we used results of Tang and Lu [36] , who derived very accurate (approximate) 
explicit analytic expressions for Kd{p) and Zd{p). Using these, one can explicitly evaluate the density derivative of 
Td{r), thus ensuring global self-consistency for the hard core reference system of the HCYF. 
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FIG. 1. Dependence of the compressibility factor PP/ p on density p* and temperature T* for the Yukawa fluid. Open 
symbols represent the results of MD simulations [30] and filled symbols are the predictions of the TIM approximation. Lines 
are guides to the eye. In the inset: dependence of the self-consistency parameter a in the TIM approach on density p* and 
temperature T* . 
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FIG. 2. Density correlation functions of the hard core Yukawa fluid. From top to bottom: (p* = 0.8, T* = 
(p* = 0.4, T* = 1.25), (p* = 0.05, r* = 1.0). MC simulations performed in this work. For the sake of clarity, 
have been shifted in the vertical direction. 
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FIG. 3. Density correlation functions of the Lennard- Jones fluid. From top to bottom: (p* = 0.54, T* = 1. 
(p* = 0.45, T* = 1.35), (p* = 0.1, T* = 1.35). MD simulations are taken from [31]. For the sake of clarity, curves 
been shifted in the vertical direction. 
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